Mesoscopic lattice Boltzmann modeling of flowing soft systems 
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A mesoscopic multi-component lattice Boltzmann model with short-range repulsion between dif- 
ferent species and short/mid-ranged attractive/repulsive interactions between like-molecules is in- 
troduced. The interplay between these composite interactions gives rise to a rich configurational 
dynamics of the density field, exhibiting many features of disordered liquid dispersions (micro- 
emulsions) and soft-glassy materials, such as long-time relaxation due to caging effects, anomalous 
enhanced viscosity, ageing effects under moderate shear and flow above a critical shear rate. 
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The rheology of flowing soft systems, such as emul- 
sions, foams, gels, slurries, colloidal glasses and related 
fluids, is a fast-growing sector of modern non-equilibrium 
thermodynamics, with many applications in material sci- 
ence, chemistry and biology [J!j. These materials exhibit a 
number of distinctive features, such as long-time relax- 
ation, anomalous viscosity, aging behaviour, whose quan- 
titative description is likely to require profound exten- 
sions of non-equilibrium statistical mechanics. The study 
of these phenomena sets a pressing challenge for com- 
puter simulation as well, since characteristic time-lenghts 
of disordered fluids can escalade tens of decades over the 
molecular time scales. To date, the most credited tech- 
niques for computational studies of these complex flowing 
materials are Molecular Dynamics and Monte Carlo sim- 
ulations [2|. Molecular dynamics in principle provides a 
fully ab-initio description of the system, but it is lim- 
ited to space-time scales significantly shorter than ex- 
perimental ones. Monte Carlo methods are less affected 
by these limitations, but they are bound to deal with 
equilibrium states. As a result, neither MD nor MC can 
easily take into account the non-equilibrium dynamics of 
complex flowing materials, such as micro-emulsions, on 
space-time scales of hydrodynamic interest. In the last 
decade, a new class of mesoscopic methods, based on min- 
imal lattice formulations of Boltzmann's kinetic equation, 
have captured significant interest as an efficient alterna- 
tive to continuum methods based on the discretization of 
the Navier-Stokes equations for non- ideal fluids [3(. To 
date, a very popular such mesoscopic technique is the so- 
called pseudo-potential-Lattice-Boltzmann (LB) method, 
developed over a decade ago by Shan and Chen (SC) [J] . 
In the SC method, potential energy interactions are rep- 
resented through a density-dependent mean-field pseudo- 
potential, *f?[p], and phase separation is achieved by im- 
posing a short-range attraction between the light and 
dense phases. In this Letter, we provide the first numeri- 
cal evidence that a suitably extended, two-species, meso- 
scopic lattice Boltzmann model is capable of reproduc- 
ing many features of soft-glassy (micro-emulsions) , such 
as structural arrest, anomalous viscosity, cage-effects and 
ageing under shear. The key feature of our model is the 
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FIG. 1: The two components A and B interact via a repul- 
sive pseudo-potential, which supports a surface tension oab- 
Moreover, each component experiences an attractive interac- 
tion in the first Brillouin zone and a repulsive one acting on 
both Brillouin zones. Each of these interactions can be tuned 
through a separate coupling constant. 



capability to investigate the rheology of these systems on 
space-time scales of hydrodynamic interest. 

The kinetic lattice Boltzmann equation takes the fol- 
lowing form [3J: 



fi.(r+8i,t+At)-Mr,t) = - — [f ls (r,t)-f^ q) (f,t)]+F is At 

where fi S is the probability of finding a particle of 
species s at site f and time t, moving along the ith 
lattice direction defined by the discrete speeds ci with 
i = 0, ...,b. The left hand-side of ([T]) stands for molec- 
ular free-streaming, whereas the right-hand side rep- 
resents the time relaxation (due to collisions) towards 
local Maxwellian equilibrium on a time scale t s and 
Fi s represents the volumetric body force due to inter- 
molecular (pseudo)-potential interactions. The pseudo- 
potential force within each species consists of an attrac- 
tive component, acting only on the first Brillouin region 
(belt), and a repulsive one acting on both belts, whereas 
the force between species is short-ranged and repulsive: 




FIG. 2: The nominal surface tension of the two-component 
fluid as a function of the reference density po (left panel). 
The coupling parameters are Ga\ = —12.55, Ga2 = 11.70, 
Gbi = -11.80, G B 2 = 10.95, G x = 0.58. By increasing the 
reference density po, the relative strength of interspecies re- 
pulsion is weakened as compared to other interactins, thus 
leading to a net decrease of the surface tension. Above a 
critical value, p cr it ~ 0.72, the nominal surface tension turns 
negative. Right panel: A typical snapshot of the density field 
Pa at time 2 x 210 6 (lattice units) at resolution 512 2 . Be- 
cause of the small but positive value of a, the system proves 
capable of supporting fairly complex metastable density con- 
figurations. 
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In the above, the indices k = 1,2 refer to the first and 
second Brillouin zones in the lattice (belts, for simplic- 
ity), c ki: pki,Wi are the corresponding discrete speeds 
and associated weights. Gab = G ss > — G s ' s , s' ^ s, is 
the cross-coupling between species, p$ a reference den- 
sity to be defined shortly and, finally rki = f+CkiAt are 
the displacements along the i-the direction in the fc-th 
belt. These interactions are sketched in Figure [TJ Note 
that positive(negative) G code for repulsion(attraction) 
respectively. Our model is reminiscent of the potentials 
used to investigate arrested phase-separation and struc- 
tural arrest in charged-colloidal systems ffj}, and also 
bears similarities to the NNN (next-to-nearest-neighbor) 
frustrated lattice spin models Q. As compared with 
lattice spin models, in our case a high lattice connec- 
tivity is required to ensure compliance with macroscopic 
non-ideal hydrodynamics, particularly the isotropy of po- 
tential energy interactions, which lies at the heart of 



the complex rheology to be discussed in this work. To 
this purpose, the first belt is discretized with 9 speeds 
(&i = 8), while the second with 16 (b 2 = 16) for a 
total of b — 24 connections. The weights are chosen 
in such a way as to fulfill the following normalization 

constraints [1]: J2iLo w i = J2iLoP^ + T,iLoP*2 = 1; 

Eiio^ c ? = E^oKi4+E^oKi4 = <2, el = 1/3 be- 
ing the lattice sound speed. The pseudo-potential ^ s (r) 
is taken in the form first suggested by Shan and Chen 
j4|, ^ s [p] = \/po(l — e~ p / po ), where po marks the density 
value at which non ideal-effects come into play. Following 
[9[ , Taylor expansion of ((2]) to fourth-order in At delivers 
the non- ideal pressure tensor P a p(r; t), namely 



P a p = [c 2 sP +-c 2 s G A i^ 2 A + G B1 ^ 2 B + ciU}S afj ~ci la p (3) 
where greek indices run over spatial dimensions and: 
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In the above equations, we have introduced the effec- 
tive couplings G a i = G" + G r s and G s2 = G s i + ^G s2 , 
s = A, B, respectively. The non ideal pressure splits 
into a local (bulk) and non-local (surface) contributions, 
which fix the surface tension er of the model. It is cru- 
cial to appreciate that the value of a can be tuned by 
changing the reference density po- The repulsive intra- 
species force F£ (proportional to y/po) acts against the 
inter-species repulsive force F x (proportional to 1/po). 
Thus, for small po, F x dominates and a complete sep- 
aration between the two fluids is expected. This is the 
case of large and positive a. On the other hand, for large 
POi cr becomes smaller and even negative. In figure [2] (left 
panel) , we show a as a function of the reference density 
P01 as obtained through a standard Laplace test on a sin- 
gle bubble configuration. As anticipated, the surface ten- 
sion decreases at increasing pa and becomes negative be- 
yond a given threshold, po > 0.71. In this work we shall 
be concerned only with the case of small and positive 
a. In the following, we shall discuss numerical simula- 
tions with random initial conditions for the two densities 
Pa and ps- After a short transient, the interfacial area 
reaches its maximum value and progressively tends to de- 
crease due to the effect of surface tension which drives the 
system towards a minimum-interface configuration. In 
the long term, this minimum-area tendency would lead 
to the complete separation between components A and 
B, with a single interface between two separate bulk com- 
ponents. However, such a tendency is frustrated (hence, 
strongly retarded) by the the complex interplay between 
repulsive (short-range inter-species and mid-range intra- 
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FIG. 3: The response function R and the surface indicator 
Iab as a function of time, (expressed in units of 10 3 LB time- 
steps). The forcing is Uo — 0.1, the domain is 128 and 
the other parameters are defined in the text. The sharp de- 
crease of the response function in the initial stage indicates 
the structural arrest of the system, associated with an anoma- 
lous enhancement of the flow viscosity, about four orders of 
magnitude above the molecular value. In the left panel, cages 
are present, which manage to "trap" micro-structures inside. 
In the right panel, the cages break down and the system is 
now able to flow again. 



species) and attractive (short-range intra-species) inter- 
actions. The final result is a rich configurational dy- 
namics of the density field, as the one shown in figure 
@ right panel. In order to investigate the rheological 
properties of the composite LB fluid, we put the sys- 
tem under a shear flow U x (x,y) = Uosin(ky), U y = 0, 
and measure the response function R = M- = Q, where 
U = J2 V U(x, y)/N y and v defines the effective viscos- 
ity. Under normal flow conditions, R = 1, so that R <C 1 
provides a direct signal of enhanced viscosity and eventu- 
ally, structural arrest. The main coupling parameters are 
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Gab = +0.58, and po = 0.7. These parameters corre- 
spond to both species in the dense phase, with no phase 
transition, hence they can be regarded as descriptive of 
glassy micro/nanoemulsions, namely a dispersion of liq- 
uid within another, immiscible liquid. By letting each 
species undergo phase transitions between a dense and 
light phase, the same model could describe foamy mate- 
rials as well. The simulations are performed mostly on 
a grid 128 2 (except the one reported in figure ^ up to 
5 x 10 6 LB time-steps. In figure[3j we show the time evo- 
lution of the response i?, as well as an indicator of the 
interface area, Iab — J2 X y ^ PA • Vps. From this figure, 
we appreciate a very dramatic drop of the flow speed in 
the initial stage of the evolution, corresponding to a very 
substantial enhancement of the fluid viscosity (about four 
orders). The system remains in this 'arrested' state for a 
very long time, over three millions timesteps, until it sud- 
denly starts to regain its initial velocity through a bu mpy 
dynamics, characterized by a series of sudden jumps [10(. 
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FIG. 4: Ageing of the system. Correlation function cor- 
responding to different waiting times t w (t w — 5 10 4 , red 
squares, t m — 2 10 5 , green circles and t m = 310 5 , blue tri- 
angles) with shear stress Uo = 0.02. In the inset, we show 
the correlation function for t w — 3 10 J and Uo — 0.03: with 
increasing shear stress the structural arrest disappears. 



These viscosity jumps signal 'plastic events', whereby the 
system manages to break the density locks (cages) which 
blocked the flow in the initial phase. As a result, the 
system progressively regains its capability to flow. These 
plastic events are also recorded by the time trace of the 
interface area Iab, which exhibits an alternate sequence 
of plateaux followed by sudden down-jumps, the latter 
being responsible for the overall reduction of the inter- 
face area as time unfolds. Visual inspection of the fluid 
morphology confirms this picture. In the top panels of 
figure [31 we show the density field in an arrested state at 
time t = 2 10 5 (top-left) and in a flowing state, t = 5 10 5 
(top-right). The left figure clearly reveals the existence of 
"cages" in the density field configuration, which entrap 
the fluid inside and consequently block its net macro- 
scopic motion. Inter-domain relaxation can only take 
place in response to 'global moves' of the density field, 
i.e. the "cage" rupture. 

Due to the mesoscopic nature of the present model, 
the rupture of a single cage in the LB simulation corre- 
sponds to a large collection of atomistic ruptures, and 
consequently it leads to observable effects in terms of 
structural arrest the system. To the best of our knowl- 
edge, this the first time that such an effect is observed 
by means of a mesoscopic lattice Boltzmann model. 

To be noted that the use of high-order lattices (24- 
speeds) is instrumental to this program, since, by se- 
curing the isotropy of lattice tensors up to 8th order, 
such lattice permits to minimize spurious effects on the 
non-ideal hydrodynamic forces acting upon the discrete 
lattice fluid |{|. We next inspect another typical phe- 
nomenon of soft-glassy matter, namely ageing. To this 
purpose, following upon the spin-glass literature |ll| . we 
define the order parameter 4> = (pa — Pb) and compute 
its overlap, defined through the autocorrelation of this 



order parameter: 
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where t w is the waiting time, r is the time lapse be- 
tween the two density configurations and brackets stand 
for averaging over an ensemble of realizations. In fig- 
ure IU we show the correlation function corresponding 
to different waiting times t w (t w — 5 10 4 , red squares, 
t w = 2 10 5 , green circles and t w = 310 5 , blue trian- 
gles) for shear stress Uo — 0.02. From this figure, ageing 
effects are clearly visible, in the form of a slower than 
exponential decay of the correlation function, which sat- 
urates to a non-zero value in the long-time limit (broken 
ergodicity). In the inset of the same figure, we show 
the correlation function for t w = 3 10 5 and Uo = 0.03: 
with increasing shear stress the structural arrest disap- 
pears, which is one of the most distinctive features of 
flowing soft-glassy materials [12] • The main advantage of 
the present lattice mesoscopic approach is to give access 
to hydrodynamic scales within a very affordable com- 
putational budget. With reference to micro-emulsions 
(say water and oil), we note that the presence of sur- 
factants usually gives rise to microscopic structures of 
the order of 50 nm in size [14j. These can be likened to 
the 'blobs' observed in our simulations. With reference 
to liquid water, we have v ~ 10~ 6 (to 2 /s), which can 
be used to obtain a physical measure of the time step 
At ~ 4 ps, about three orders of magnitude larger than 
the typical timestep used in Molecular Dynamics. As a 



result, a five-million time-step LB simulation spans about 
20 microseconds in physical time. Since the present LB 
method is easily amenable to parallel computing, paral- 
lel implementations will permit to track the time evo- 
lution of three-dimensional micro-emulsions of tens of 
microns in size, over time spans close to the millisec- 
ond, i.e. at space-time scales of hydrodynamic relevance. 
Summarizing, we have provided the first numerical ev- 
idence that a two-species mesoscopic lattice Boltzmann 
model with mid-range repulsion between like-molecules 
and short-range repulsion between different ones, is ca- 
pable of reproducing many distinctive features of soft ma- 
terial behaviour, such as slow-relaxation, anomalous en- 
hanced viscosity, caging effects and aging under shear. 
The present lattice kinetic model caters for this very rich 
physical picture at a computational cost only marginally 
exceeding the one for a simple fluid. As a result, it is 
hoped that it can be used as an alternative/complement 
to MonteCarlo and/or Molecular Dynamics, for future 
investigations of the non- equilibrium rheology of a broad 
class of flowing disordered materials, such as microemul- 
sions, foams and slurries, on space and time scales of 
experimental interest. 
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